Astronomy &: Astrophysics manuscript no. CelesteBlazarAA 


February 5, 2008 


(DOI: will be inserted by hand later) 





Mrk 421, Mrk 501, and lES 1426+428 at 100 GeV with the 

CELESTE Cherenkov Telescope 

D. A. Smith^, E. Brion^*, R. Britto^, P. Bruel^, J. Bussons Gordo^, D. Dumora\ E. Durand\ 
P. Eschstruth^, P. Espigat^, J. Holder'^, A. Jacholkowska^ , J. Lavalle^, R. Le Gallou^, B. Lott^, 
H. Manseri^, F. Miinz^, E. Nuss^, F. Piron^, R. C. Rannot^'*^, T. Reposeur^ and T. Sako^ 

^ Centre d'etudes nucleaires de Bordeaux Gradignan - CENBG UMR 5797 CNRS/IN2P3 - Universite Bordeaux 1, 

Chemin du Solarium - BP120 33175 Gradignan cedex, France 
^ Laboratoire de physique theorique et astroparticules, Universite Montpellier 2, 34095, France (UMR 5207 

CNRS/IN2P3) 

^ Laboratoire Leprince-Ringuet, Ecole Polytechnique, 91128 Palaiseau, France (UMR 7638 CNRS/IN2P3) 
" Laboratoire de I'accelerateur lineaire, Universite Paris 11, 91898 Orsay, France (UMR 8607 CNRS/IN2P3) 
^ Laboratoire d'astroparticule et cosmologie. College de France, 75231 Paris, France (UMR 7164 CNRS/IN2P3) 
^ Astrophysical Sciences Division, Bhabha Atomic Research Centre, Mumbai-400 085, India 

version: 19 August 2006 

Abstract. We have measured the gamma-ray fluxes of the blazars Mrk 421 and Mrk 501 in the energy range 
between 50 and 350 GeV (1.2 to 8.3 x lO^'^ Hz). The detector, called CELESTE, used first 40, then 53 hehostats 
of the former solar facility "Themis" in the French Pyrenees to collect Cherenkov light generated in atmospheric 
particle cascades. The signal from Mrk 421 is often strong. We compare its flux with previously published multi- 
wavelength studies and infer that we are straddling the high energy peak of the spectral energy distribution. The 
signal from Mrk 501 in 2000 was weak (3.4 a). We obtain an upper limit on the flux from lES 1426-1-428 of less 
than half that of the Crab flux near 100 GeV. The data analysis and understanding of systematic biases have 
improved compared to previous work, increasing the detector's sensitivity. 
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1. Introduction 

' Measurements of high energy emission from active galac- 
. tic nuclei (AGN) give insights into a variety of open 
' problems, such as the nature of the AGNs themselves, 
and the density and evolution of the extragalactic dif- 
fuse infrared background. After GLAST is launched 
in 2007 IjGehrels fc Michelson 1999|l . the large number 
of high galactic latitude GeV gamma-ray sources to 
be seen will make AGNs become a background for 
searches for new classes of emitters - insufficient un- 
derstanding of high energy AGNs might just limit po- 
tential glimpses of a hidden Universe. Several AGNs 
of the "blazar" class emit above 250 GeV and have 
been detected by atmospheric Cherenkov imaging tele- 
scopes. They are called HBLs, for high energy BL 
Lac objects. Of th ese, Mrk 421 (fbtazejowski et al. 2005| 
IPiron et al. 2001al lAharonian et al. 2003|l and Mrk 501 
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l|Dwek fc Krennrich 2005L |Djannati-Atai et al. 1997| ) are 

the brightest . 

The Synchrotron Self Compton model ("SSC") is 
consistent with much of the multiwavelength blazar 
data. For discussion of these models as applied to 
the blazars discussed in this paper, see for example 
UTavecchio et al. 1998)1 or ( |Katarzynski et al. 200 1| ). In 
its simplest variant, relativistic electrons generate syn- 
chrotron photons with a vFi, peak in the optical to X- 
ray range, and then up-scatter these same photons to 
produce the peak seen in the GeV range. The inverse 
Compton peak position can constrain the physical param- 
eters of the model. Many, if not most, observed charac- 
teristics of blazars remain unexplained. Examples are the 
relation between variations at different wavelengths, as 
highlighted by the "orphan" gamma ray flares observed 
for lES1959-(-650 UKrawczynski et al. 2004| ), and the role 
that proton acceleration might play. Progress requires si- 
multaneous measurements across the spectrum. 

The CELESTE telescope was one of the first 
instruments sensitive near the Compton peak 
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for these blazars, beyond the energy range of 
the relatively large sample of EGRET blazars 
( |Mukherjee et al. 1997| IHartman et al. 1999|l . Mrk 421 
was studied with EGRET, but Mrk 501 was vis- 
ible only during a flare IjKataoka et al. 1999|l . 
while lES 1426+428 was undetected by EGRET, 
although well-measured at TeV energies (e.g. 
|Djannati-Atai et al. 2002) IHoran et al. 2002|l . A key 
motivation for our experiment was to provide blazar 
measurements in an unexplored energy window. 

This paper presents the first observations of three 
blazars below 100 GeV, by a solar heliostat array, and is 
structured as follows: we describe the detector, the analy- 
sis method, and the data samples. We emphasize improve- 
ments in our knowledge of the systematic uncertainties. 
We present light curves for the three blazars as well as 
the integrated fluxes, and place the results in the context 
of SSC models. 

2. The CELESTE gamma-ray Telescope 

The CELESTE detector is located in the eastern French 
Pyrenees (N. 42.50 °, E. 1.97°, altitude 1650 m) and de- 
scribed in fGicbcIs ct al. 1998, P are et al. 2002 *1 and in 
l|de Naurois et al. 2002,1 . FigureHiUustrates the principle. 
CELESTE was dismantled in June, 2004. 

CELESTE detecte d the Crab n ebula in 1998 using a 
solar heliostat array llSmith 199811 . STACEE, very sim- 
ilar to CELEST E, reported the flux above 190 GeV 
l|Oser et al. 2001|l , followed by the flux above 60 GeV 
from CELESTE l|de Naurois et al. 2002|l . STACEE re- 
cently measured the integral flux above 140 GeV from 
Mrk 421 IjBoone et al. 2002|l and, assuming power law 
spectral indices extrapolated from higher energies, brack- 
eted the 50 to 700 GcV differential spectrum. A re- 
view of solar tower gamma-ray telescopes can be found 
in l|Smith 2005|l . Here we recall key features of the 
CELESTE detector. 

CELESTE used 53 hehostats (40 until 2001) of the 
Themis former solar facility. The CAT Cherenkov imager 
ran on the same site l|Piron et al. 200 lal and references 
therein), as did the Themistocle HJjaillon et al. 1 993) 
and ASGAT IjCoret et al. 1993|l Cherenkov sampling ar- 
rays. Each heliostat has a 54 m^ mirror on an alt-azimuth 
mount. Light from the heliostats was reflected to a sec- 
ondary optical system at the top of the 100 m tower, where 
one photomultiplier assembly (PMT) viewed each helio- 
stat. Heliostat control and data acquisition systems were 
also high in the tower. 

Data acquisition was triggered when the summed PMT 
pulse height for at least M of the N groups of heliostats 
exceeded 4.5 photoelectrons ("p.e.") per heliostat. For the 
40 heliostat array, Af = 3 of iV = 5 groups with 8 he- 
liostats each was used most often. For the 53 heliostat 
array there were N = 6 trigger groups of between 6 and 
8 heliostats each, and the trigger multiplicity was set at 
M = 3,4, or 5 for different observations. The deadtime 
in the electronic delays depends on the rates, which de- 




Fig. 1. Experimental principle: Heliostats tracking a 
source reflect Cherenkov light from atmospheric particle 
cascades to the secondary optics, photomultipliers, and 
acquisition electronics located high in the 100 meter tall 
tower. 

pend on the threshold settings, which are determined in 
part by the choice of M. We took pains to remove these 
biases UManseri 2004|l . in particular, with an offline group 
multiplicity cut based on the scaler rates. The acceptance 
shown in figure has an offline multiplicity cut of 4 of 5, 
for illustration. 

Each PMT signal was digitized at 1.06 nanoseconds 
per sample, with about 3 digital counts (dc) per photoelec- 
tron. 100 samples centered around the nominal Cherenkov 
pulse arrival time were recorded, providing event-by-event 
pedestal information that we used to determine the night 
sky background light for each channel. A second 28 sample 
window contained a timing reference pulse. Scaler rates, 
PMT anode currents, a GPS time stamp, and some mete- 
orological information rounded out the data record. 

3. Monte Carlo Simulations 

celeste's first publication, on the flux from the Crab 
nebula quoted a ~ 30 GeV threshold at the trigger level, 
and 60 GeV after analysis Ijde Naurois et al. 2002,1 . We 
estimated the uncertainty on the energy scale to be less 
than 30 %, dominated by disagreement between predic- 
tions of Monte Carlo atmospheric cascade generators, and 
disagreement between estimated and observed PMT illu- 
minations induced by bright stars. Since then we have cor- 
rected some small errors in the KASCADE Monte Carlo 
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IjKertzman &: Sembroski 1994|l and the agreement is im- 
proved (Jemc = 5 to 10 %). T he CQRSIKA Mon te Carlo 
atmospheric shower generator IjHeck et ai. 1998)l best re- 
produces our experimental observables, and we use it for 
the results presented here. 

Cherenkov photons from the shower generator are fed 
into a detailed simulation of our optics and electronics, to 
develop our data analysis methods and cuts, as well as to 
calculate the energy dependent effective area A{E). Here 
we describe improvements made to the detector simula- 
tion since l)de Na urois ct al. 2002). The acceptances will 
be discussed after the analysis cuts have been explained. 
Simulations indicate that the trigger rate is dominated by 
18 Hz of protons, with 4 Hz of helium nuclei. Photons from 
the Crab nebula trigger the detector at about 2% of the 
trigger rate. 

3.1. Optical Throughput 

Previously, the phototube illuminations when viewing 
stars were lower than predicted by the simulation. We 
reviewed the ray tracing programs, and found two fac- 
tors that dominated the discrepancy. First, the combined 
heliostat and secondary mirror reflectivities are about 
20 % worse than earlier measurements indicated ^ . Second, 
the simulated heliostat focussing was sharper than in re- 
ality. We compared the image sizes obtained from star 
"scans" used for heliostat alignment with simulated sizes, 
and modified the simulation to improve agreement. (In 
scans, we record PMT anode currents as the heliostats 
are pointed successively at a grid of points around the 
star position: see figure 7 in lPare et a,l."2nf)2|l . 

Figure |2] illustrates a consequence of the ray tracing 
modifications HBrion 2004|l . When tracking Mrk 421, the 
star 51 UMa (HD 95934, blue magnitude Mb = 6.16) 
falls within the ±5 nirad optical field-of-view for heliostats 
near the center of the array. The On-sourcc PMT anode 
currents are thus larger (ion ~ 13 M-A) than the Off- 
source currents which are due only to the diffuse night 
sky light («off ~ 10 M-A). Using the PMT gains g, the 
current differences yield the star-induced illuminations 
b = («On — *Off)/ffe, in photoelectrons per second (e is 
the electron charge). The figure shows that measured illu- 
minations vary with the pointing direction, but less than 
was predicted by the original simulation. The simulations 
include the stars in the On and Off fields that contribute 
at least 5 % as much as 51 UMa, that is, Mb < 9.4. We 
remark that the gamma-ray energy scale predicted by the 
Monte Carlo simulation shifts less than figure |21 suggests, 
since an atmospheric shower is a light source of angular 
extent comparable to the detector field-of-view, and so the 
effect of heliostat "defocussing" is less violent for showers 
than it is for the point-like stars. 



^ We thank Professor Ladislav Rob of Charles University, 
Prague, and his colleagues at Compas Ltd., Turnov, Czech 
Republic, for the recent measurements. 



A few observations constrain the uncertainty in the 
optical throughput. Photometry using star scans yields 
20 % channel-to-channel dispersion compared to the sim- 
ulated values. A census of damaged heliostat mirrors, and 
the phototube manufacturer's measurements of the pho- 
tocathode responses, combined together account for half 
of the dispersion. We attribute the rest to errors in the 
gains, and focussing anomalies for individual heliostats. 
However, what matters most is the sums over the 8 he- 
liostats of the trigger groups: the total variation over the 
5 groups is ±15 %, and the rms is less than 10 %. We 
conclude that the overall effect on the gamma-ray energy 
scale of the uncertainty in the optical throughput Seopt is 
less than 10 %. 
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Fig. 2. Example of phototube illuminations induced by 
the star 51 UMa in the field of view of Mrk 421, ver- 
sus hour angle, for hehostat C07. Squares & crosses: For 
the 40 pairs of On- and Off-source "double-pointing" data 
runs passing current and rate stability criteria, the quan- 
tity shown is the difference of the average current dur- 
ing each run divided by the PMT gain g, (ion — ios)/g- 
For the crosses, the trigger rates were stable but lower 
than usual: optical transmission was low, presumably 
due to clouds or frost. Circles: Original detector simu- 
lation. Stars: Simulation results, after "defocussing" the 
heliostats and degrading the mirror reflectivities as de- 
scribed in the text. The curves are to guide the eye. 

Since Ijde Naurois et al. 2002|l . the electronics simula- 
tion has been completely re-written, leading to improved 
agreement between predicted and observed signals, in the 
trigger, and for single channels. The good agreement for 
the final analysis variable, described below, is the result 
of the care taken to model the elements of the complex 
electronics chain. 

3.2. Atmospheric Transmission 

The amount of Cherenkov light that reaches the 
ground varies between different sites and seasons de- 
pending mainly on the aerosol content versus altitude 
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Source SP DP SPV All pairs < 15 Hz On hours 



Mrk 421 


61 


44 


28 


133 


20 


% 


38.9 


Mrk 501 


10 


25 




35 


50 


% 


10.3 


lES 1426+428 






33 


33 


40 


% 


8.8 



Table 1. Number of data pairs, after applying stability criteria to the PMT anode currents and to the trigger rates. 
Single Pointing (SP), Doubling Pointing (DP), and Single Pointing with Veto (SPV) are explained in the text. The 
fraction of the data set with a cosmic ray trigger rate below 15 Hz, which we ascribe to seasonal increases in atmospheric 
extinction, is also listed. 



HBernlohr 2000|l . The same also holds for atmospheric ex- 
tinction in stellar photometry ( |Hayes fc Latham 197 5). 
The effect on CELESTE is particularly striking: the 
25 Hz cosmic ray trigger rates for winter sources (e.g. Crab 
and Mrk 421) decrease to under 15 Hz every year in spring, 
affecting sources like Mrk 501 and lES 1426+428. Table □ 
lists the fraction of the data sets above and below 15 Hz. 

The CORSIKA Monte Carlo provides a choice of 
tabulated wavelength and altitude dependent extinction 
curves. The default contains some aerosols. The Palomar 
curve from (Ha yes fc Latham 1975| ) used for our stellar 
photometry studies contains fewer aerosols, for an inte- 
grated transmission above Themis at 400 nm 6 % greater 
than for the default CORSIKA curve. 

An amateur telescope with filters and a CCD cam- 
era was used to measure the total integrated atmospheric 
transmission during some nights when the trigger rates 
were above 20 Hz, using the apparent brightness of sev- 
eral stars as a function of their zenith angles (Bouguer 
method), at three wavelengths (blue: 390-490 nm, green: 
500-560 nm, and red: 610-660 nm). The Cherenkov light 
recorded by CELESTE is mainly in the blue range. The 
CCD results show that the sky at Themis corresponds to 
the Palomar extinction curve. We made a CORSIKA ex- 
tinction table corresponding to the Palomar curve which 
we use to simulate good running conditions. 

To calculate the acceptances for data acquired with 
trigger rates < 15 Hz we made another CORSIKA extinc- 
tion table that we call "dusty" . We increased the aerosol 
contribution to the Palomar curve so that the integrated 
transmission above Themis at 400 nm is 60 % less. This 
choice comes from two observations. The first appears in 
Figure |2 where the phototube illumination due to the 
bright star near Mrk 421 in good conditions is halved for 
data acquired with a low trigger rate. The second is that 
proton shower simulations using the nominal and "dusty" 
extinction tables reproduced the experimental nominal 
and degraded trigger rates, respectively. Figure [SJ dis- 
cussed below, shows the acceptance curves obtained for 
the two atmospheres, used to convert our gamma-ray rates 
to fluxes. 

4. Data Samples 

Maximum Cherenkov emission for 100 GeV 7-ray show- 
ers occurs around 11 km above the site. For CELESTE'S 
oldest data, the heliostats were aimed at this height in 



the direction of the source, to maximize Cherenkov light 
collection. This is called Single Pointing data (SP). 

SP provides the lowest energy threshold, but also de- 
creases the detector's sensitive area. For one season we 
took data in a configuration called Double Pointing (DP), 
where half the heliostats point at 11 km and the rest at 
25 km. This increases the area and improves cosmic ray 
rejection, while raising the energy threshold only slightly. 

Our ±5 mrad field-of-view matches the apparent size of 
a gamma shower, to optimize the ratio of Cherenkov light 
to night sky background light. Tails of cosmic ray showers 
are outside the field-of-view, reducing background at the 
trigger level but weakening offline rejection. In 2001 we 
added 13 heliostats to CELESTE, and thereafter pointed 
41 heliostats at the 11 km point, and the remaining 12 
at a ring 150 meters around the point. The 12 "veto" 
heliostats can see shower tails outside the central field-of- 
view, enhancing cosmic ray rejection IjManseri 2004|l . We 
call this Single Pointing with Veto data (SPV). 

Ultimately the performance of the analysis was such 
that the veto rejection was not used. But the 12 veto he- 
liostats around the edge of the heliostat field led us to 
restructure the trigger into M ~ 6 groups. The changed 
trigger topology affects the shape of the recorded showers, 
leading to different energy dependent acceptances. 

Source observations ( "On" ) last about 20 minutes and 
are followed or preceded by an observation at the same 
declination, offset by 20 minutes in right ascension. "Off" 
data gives a reference for the cosmic ray background - the 
signal is the difference between On and Off, after analy- 
sis cuts. However, changes in sky conditions between On 
and Off can cause differences having nothing to do with 
gamma-rays: we select stable On-Off pairs with a cut on 
the value obtained from a least-squares fit to a con- 
stant PMT anode current versus time, and to the trigger 
rate versus time, for both runs. We also rejected pairs with 
very low data rates (< 10 Hz), or anomalously high trig- 
ger dead times. The weather at Themis is fickle, and we 
rejected over half of our data. 

Table^lists the data sets, recorded during clear, moon- 
less nights from December 1999 through May 2004, after 
data selection. 



CELESTE: Mrk 421, Mrk 501, and lES 1426+428 at 100 GeV 



5 




Fig. 3. LEFT: Summed Cherenkov pulse - a 100 ns Flash ADC window, centered at the time an ideal spherical 
wavefront 11 km directly above the site would reach the photomultiplier tubes after reflection from the heliostats, 
is recorded for each channel. In analysis, we sum the 40 (SP and DP data) or 41 (SPV data) channels repeatedly 
over a grid of hypothetical shower core positions. RIGHT: The H/W values for each position of the grid, for a single 
simulated event. The grid position with (-ff/W^)max is a good (±15 m) estimator of where the gamma-ray would have 
hit the ground. The analysis variable ^ — {H/W)2()o m/(^^/W^)inax uses {H/W)2oo m, which is the value averaged over 
a ring 200 meters from the grid center. The vertical axis is in ADC counts per nanosecond. 



5. Data Analysis 

5.1. Event reconstruction and background rejection 

Since Ijde Naurois et al. 2002|l we have changed our anal- 
ysis strategy, as described below. The Flash ADC data 
quality improved, as did the channel-to-channel timing off- 
sets, which both help wavefront reconstruction. Our pro- 
cedures to remove the biases induced by night sky back- 
ground differences at the trigger level ( "software trigger" ) 
and in the FADC data ("padding") have been refined. 
Overall sensitivity has more than doubled. 

Just above the trigger threshold, the Cherenkov sig- 
nal in one channel is comparable to the fluctuations of 
the night sky background light. Summing the channels 
greatly improves the signal-to-noise ratio, after correcting 
for the propagation delay to each heliostat, assuming a 
spherical Cherenkov wavefront. Lacking knowledge of the 
shower core position, we scan the plane at 11 km and 
evaluate the ratio H/W for each point in a grid, where 
H and W are the height and the width of the summed 
signal (see Figure OJ. The grid position giving the largest 
value of the ratio H/W is a measure of the shower core 
position. The wavefront sphericity correlates with how 
much H/W decreases 200 m away from the shower core. 
We measure the relative decrease with a variable called 
^ = {H/W)2oo m/{H/W)ma.7,, showu in Figure H for an 
Off observation, Mrk 421 On-Off data, and for a simula- 
tion of a 7-ray spectrum. The Cherenkov wavefront from 
7-ray induced showers is, on average, more spherical than 
for showers initiated by charged cosmic rays and, as ex- 
pected, 7-ray showers have smaller ^ than showers initi- 
ated by charged cosmic rays. The method and prelimi- 



nary results were shown in IjBruel 2004|l and detailed in 
IjManseri 2004|l . 

The ^ distributions change for different heliostat point- 
ings. This article concerns three blazars with similar de- 
clinations, all transiting near zenith at Themis, so the only 
changes are the pointing configuration and the hour an- 
gle. We studied simulations and the strong signals from 
Mrk 421 and the Crab to explore the ^ cuts l|Brion 2005|l . 
Statistical significance (sensitivity) is optimized by cutting 
a bit above the peak in the gamma ray ^ distribution, e.g. 
^ < 0.35 in Figure 01 For other pointings, the peak lies 
between ^ ~ 0.25 and ~ 0.35. The acceptance depends 
in part on the efficiency of the ^ cut, obtained from the 
Monte Carlo. This efficiency is related to the integral of 
the ^ distribution, a steep function of ^ near the peak. 
Thus we choose to cut beyond the peak, at ^ < 0.4, for all 
data sets, which increases the statistical uncertainty, but 
decreases the acceptance uncertainty Seacc to less than 
10%, and is simpler than having an array of cut values for 
the various pointings. 

5.2. Acceptances and Performance 

For a gamma-ray source with differential spectrum (j^iE) 
the number N of detected gammas will be 

/>oo 

N = T A{E)(j){E)dE 
Jo 

where T is the On-source time and A{E) is the en- 
ergy dependent gamma-ray collection area obtained from 
the Monte Carlo calculations. We have calculated A{E) 
for each heliostat configuration, and for hour angles of 
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Fig. 4. LEFT: ^ distributions for Mrk 421 for March 2004 data. Black dots show On-Off data. The broad histogram is 
Off data (i.e. cosmic ray background). The narrow histogram is simulated gamma-rays. The simulation and Off data 
are normalized to the number of On-Off events. RIGHT: Statistical significance of the excess as a function of the ^ 
cut value. 



(0, 1.0, 1.5, 2.0) hours, for each of the two atmospheres de- 
scribed above. Figure 13 shows A{E) at the trigger level, 
and after cuts, for the direction of the blazar Mrk 421 
{S = 38°), one hour after transit. (The few degree discrep- 
ancy with the other two blazar's declinations induces a 
negligible error) . For the blazar spectra we assume a power 
law, ^(i?) = kE~°'. celeste's energy range is near the 
SED high energy peak for the three blazars, where a = 2. 
We vary a between 1.8 and 2.2 to estimate the uncertainty 
due to this choice. The energy reconstruction method in 
HPiron et al. 2008(l was used to search for a signal in en- 



ergy bands in (I Lavalle et al. 2006)l but was not applied 
here. 

The combined uncertainties (Scmc ® ^^opt) < 20% af- 
fect the energy scale. Furthermore, along with Seacc they 
lead to a flux error. A remaining uncertainty comes from 
the varying atmospheric transmission. We make a first- 
order correction by grouping the data into subsets with 
trigger rates of < 15 and > 15 Hz (see Table^ for use with 
the nominal or "dusty" acceptances, A{E), as in Figure 
The maximum counting rate after cuts is nearly halved 
with the "dusty" atmosphere. The observed extrema of 
the trigger rates are 12 and 24 Hz, but the distribution 
clusters around 14 Hz and 22 Hz. We conclude that the 
overall systematic flux uncertainty is < 25%. 

To summarize CELESTE'S performance: for one hour 
of On-source Crab data near transit, analysis using ^ < 0.4 
gives 3.7 ± 0.2 gamma/min over the whole data set. The 
significance is 3.3 a for SP and DP data, and 5 <t for 
SPV. The higher sensitivity for the more recent data is 
due mainly to having six trigger groups instead of five, 
at the cost of a somewhat higher minimum energy: Table 
121 shows that the SP raw data rate is 15% higher than 
for SPV, but that analysis cuts reject more background 
for SPV data (90%) than for SP data (83%). The higher 
sensitivity also stems in part from improvements in FADC 
quality, timing adjustments, and gain settings over the 



years. Optimized cuts increase the significances by 20 %, 
and more for Mrk 421 IjErion 2005|l but are not used in 
this paper. 

Folding our energy dependent acceptance A{E) (simi- 
lar to Figure 01 with the Crab inverse Compton spectrum 
parametrized in Table 6 of (|Aharonian et al. 2004|l yields 
a maximum at 90 GeV. The rate falls to 20 % of the 
maximum below 50 GeV and above 350 GeV. Figure [HI 
shows the acceptance-corrected pair-by-pair Crab fluxes, 
as a function of date and hour angle. The stability is ade- 
quate to search for blazar variability. The power flux at 90 
GeV is vF^ = 0.89±0. 05 x 10~^° erg cm"^ s~\ higher than 
our previous result in Ijde Naurois et al. 2002|l but within 
the uncertainties, and in agreement with the spectrum in 
HAharonian et al. 2004|l . 

6. Results and Discussion 

Tables 1201 and |S1 summarize the data analysis results for 
the three blazars Mrk 421, Mrk 501, and lES 1426+428, 
respectively. Figures |H| and show the acceptance- 
corrected light curves for the blazar Mrk 421 during the 
years 2000, 2001, and 2004 seasons. Figures [13 and ITl 
show the light curves for Mrk 501 and lES 1426+428, 
respectively. The curves assume power laws with spectral 
index a = 2. There is an often strong signal from Mrk 421, 
and no signal for lES 1426+428. The excess in the Mrk 501 
data is discussed below. 



6.1. Mrk 421 

Understanding blazars requires observations at different 
wavelengths, preferably simultaneous given their vari- 
ability. Happily, much of our data coincides with pre- 
viously published results. Table |3 lists the CELESTE, 
TeV, and X-ray fluxes for the epochs with CELESTE 
data, celeste's richest data taking period was March 
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Fig. 5. LEFT: Energy dependent gamma-ray collection area A{E), at the trigger level (round markers), and after 
analysis cuts (triangles), for sources culminating near zenith, 1 hour after transit, for 11 km single pointing (SP). The 
curves are splines. The two solid curves used the low-aerosol atmosphere while the two dashed curves used the "dusty" 
atmosphere (see text). The analysis trigger group multiplicity cut for these plots was M = 4 of the N = 5 groups. 
RIGHT: A{E) multiplied by a typical flux spectrum, to give the relative counting rate A{E)/E'^ . Integral fluxes 
are quoted at the energy corresponding to the maximum rate, while the energy range of 50 to 350 GeV corresponds 
to 20% of the maximum rate. 




51400 51600 51800 52000 52200 52400 52600 52800 53000 53200 -2-10123 

(MJD) Hour Angle(hours) 



Fig. 6. LEFT: Crab fluxes for single data runs over the 5 years that CELESTE ran. SP pointing was mainly used 
during the first year, DP for the 2nd year, and SPV for the last. The integral flux above 90 GeV is 0.62 ± 0.04 x 10~^ 
photons/cm^/s averaged over the three good Crab seasons. RIGHT: Same, as a function of the hour angle of the data 
runs. After acceptance corrections the observed flux is steady. 



Number of events Significance Signal-to-noise Excess 



Data set 


Cut 




A^On 


iVoff 


iVon - iVofI 




[%] [7/min] 


SP (17.8 h) 


Raw Data 
Software trig 
All cuts 


ger 


1 575 396 
1 056 705 
276 057 


1557192 
1044 811 
270113 


18 204 
11894 
5 944 


7.1 
7.1 


1.1 11.1 

2.2 5.6 


DP (12.9 h) 


Raw Data 
Software trig 
All cuts 


ger 


983 550 
532 883 
113 507 


961 257 
523 661 
106 904 


22 293 
9 222 
6 603 


8.1 
12.8 


1.8 11.9 
6.2 8.5 


SPV (8.2 h) 


Raw Data 
Software trig 
All cuts 


ger 


621 243 
410 429 
62 914 


604 926 
401 791 
57551 


16 317 
8 638 
5 363 


8.5 
13.9 


2.1 18.3 
9.3 11.8 


All data (38.9 h) 


Raw data 
Software trig 
All cuts 


ger 


3180189 
2 000 017 
452 478 


3123 375 
1 970 264 
434 568 


56 814 
29 754 
17910 


13.1 
16.9 


1.5 12.7 
4.1 7.7 



Table 2. Number of events before and after cuts, for the Mrk 421 data sets. N„ is calculated after correcting for the 
~ 80% data acquisition efficiency. 
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Fig. 7. Mrk 421 daily averages in early 2000. The CELESTE data runs from MJD 5 1572 to 51613 ( January 29 to March 
10), while epoch (1) is defined in Table|2|and discussed in the text. CAT data is from IjPiron et al. 2001a. left-hand scale) 
while the HEGRA points are from (jKrawczynsk^ t al. 20011 right-hand scale). RXTE ASM Ijsee ASM in references! 
left-hand scale) and PCA ( |KrawczyiisIa'enLr'200T[ right-hand scale) results are also shown. 



2000, with 11 of 14 consecutive nights (epoch (1), 
from 2000 February 27 to March 12, MJD 51601.96 to 
51615.09). Figure H sho ws the daily averaged data dur- 
ing this period. CAT (Piron et al. 2001a|l and HEGRA 
IjAharonian et al. 2002(1 have data for most of the same 
nights, and the four experiments with data at this epoch 
are all at a low level. The CELESTE flux averaged over 
the 2003 season, epoch (3), from February through April 
(MJD 52667.12 to 52762.94) was at an intermediate level, 
as were the other experiments. Overall, the flux variations 
observed with CELESTE track the other measurements. 

Figure ^ shows the SED for Mrk 421 taken from 
( iBTazejowski et al. 2005| ) to which we have added our re- 
sult for epoch (2) when TeV data from Cherenkov im- 
agers are available, 4 nights in February 2001 (MJD 
51956.98 to 51963.11). ASM indicates that the blazar was 
at the low end of the "medium" X-ray state as defined in 
HBtazejowski et al. 2005|). We have plotted the CAT spec- 
trum for the vear 2000 l|Piron et al. 2001a|l . which has the 
same integral flux as measured on those four nights. 

We now look more closely at the light curves. 
The detailed multi-wavelength study reported in 



HBlazejowski et al. 2005| ) makes February-March 2004 
particularly interesting (epoch 4). CELESTE measured a 
steady high state of 1.04 ± 0.1 x 10~^ cm~^s~-^ for the 
three nights of 2004 February 15 to 18 (MJD 53050.05 to 
53053.08, see Tableland Figure 0. VERITAS shows a 
steady low state, in seeming contrast to CELESTE. But 
Figure 1 of ( ^Blazejowski et al. 2005 ) shows an increasing 
optical flux on those dates, while PCA and ASM on 
RXTE show a decreasing X-ray state. This could be a 
shift of the low energy peak of the SED to lower energy. 
If the high energy SED peak also shifts during this epoch 
it could explain that CELESTE sees a high rate when 
VERITAS does not: VERITAS is above the SED peak 
while the CELESTE range straddles it. 

Similarly, and also in Figure |21 PCA and ASM show 
a declining X-ray state followed by a small rebound for 
2004 March 16 to 18 (epoch 5). For two nights of data 
(MJD 53079.91 to 53082.02) CELESTE has an average 
state of 1.4±0.2 x 10"^ cm^^s-i, but here VERITAS sees 
the decline of the flare before returning to the low state. 
CELESTE'S stability during VERITAS' changes is again 
consistent with Figure [Till that is, with an SED peak at 
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Epoch 


ASM ( 


a) 


PCA 


CAT (d) 


HEGRA (e) 


VERITAS (f) 


CELESTE (g) 


2000 Feb 4 


1.25 


± 


0.10 




189.7 ± 21.2 


11.24 ± 3.41 








0.97 


± 0.17 


(51577.99 to 51578.18) 
























2000 Feb 10 


1.51 


± 


0.21 


1.41 ±0.02 (b) 


85.3 ± 19.7 


2.43 ± 1.93 








0.48 


± 0.14 


(51583.99 to 51584.18) 
























(1) 2000 Feb 27 to March 12 


0.39 


± 


0.06 




48.3 ± 10.9 


3.9 ± 0.8 








0.29 


±0.07 


(51601.96 to 51615.09) 
























(2) 2001 Feb (4 nights) 


2.98 


± 


1.57 




208.7 ± 10.5 


22.5 ± 0.8 


4.35 


± 


0.29 


0.93 


± 0.09 


(51956.98 to 51963.11) 
























2001 April 14 to 17 


1.26 


± 


0.17 




131.7 ± 16.3 


18.7 ±3.1 


4.65 


± 


0.22 


0.34 


±0.12 


(52012.90 to 52015.97) 
























(3) 2003 Feb to April 


0.70 


± 


0.03 


17.3 ± 0.02 (c) 






0.94 


± 


0.04 


0.60 


±0.12 


(52667.12 to 52762.94) 
























(4) 2004 Feb 15 to 18 
(53050.05 to 53053.08) 

(5) 2004 March 16 to 18 
(53079.91 to 53082.02) 


2.47 ±0.23 
1.74 ±0.31 


39.3 ±0.1 (c) 

48.4 ±0.2 (c) 






1.66 ±0.20 
3.37 ±0.17 


1.04 ±0.10 
0.71 ±0.09 



Table 3. Multiwavelength fluxes from Mrk 421 coincident with CELESTE data. The epoch numbers (1) through (5) are 
referred to in the text and figures, (a) ASM counts per second Ijsee ASM in referencesjl . RXTE PCA 2-10 keV fluxes are 
in (b) 10~^° ergcm"^ s~^, 3-20 keV ( |Krawczynski et al. 2001| ), and in (c) counts per second ( |Blazejowski et al. 2005| ). 
(d) CAT fluxes are in 10" photons cm "^ s~\ E > 250 GeV (Piron ct al. 2001al 2 001b). (e) HEGRA flu xes are in 
10-12 photons cm-2s-\ E > 1000 GeV ( |Krawczynski et al. 2001, for the 2000 data. lAharonian et al. 200 31 for 2001). 
(f) VERITAS fluxes above 300 GeV are in ga mma/min, where the Crab gave 2.40 (2.93) gamma/min for 02/03 (03/04) 
l|Holder et al. 20011 [Biazejowski et al. 2005| ). (g) CELESTE fluxes are in 10"^ photons cm^^ g-i, E > 100 GcV. 



100 GeV. Recent spectral work on Mrk 421 by STACEE 
also seems to indicate that our 50 to 350 GeV energy range 
matches a high energy SED peak IjCarson 2005|l . 

6.2. Mrk 501 

The statistical significance of the excess for all the Mrk 501 
data is weak (2.9 a), as shown in Table 01 The same 
^ < 0.4 cut as for Mrk 421 was applied. The light curve 
in Figure^] shows that the excess comes mostly from the 
year 2000, with 3.4 a for those four months. Figure 
shows that the excess in 2000 ressembles a gamma-ray 
signal. However, changing the ^ cut values to those that 
optimize the statistical significance of the Mrk 421 sig- 
nal do not enhance this excess, suggesting that taking the 
excess to be a gamma-ray signal is delicate. 

CELESTE and the BeppoSAX X-ray satellite have 
data in common for four nights, from 2000 March 5 to 
12, when Bcppo measured between (0.76 ± 0.08) and 
(1.18 ±0.12) X 10-1° ergcm-^s-i (|Massaro et al. 2004jl . 
The ASM rate during this period is a roughly con- 
stant 0.4 ±0.1. (ASM rates are unreliable at such low 
fluxes.) Averaging the ASM count rates for the nights 
with CELESTE data gives 0.44 ± 0.06 and 0.37 ± 0.09 for 
2000 and 2001, respectively. The dispersion of the points 
is 0.5 cts/s for both seasons. Averaging ASM over longer 
periods gives essentially the same result, so that no sig- 
nificant difference in X-ray activity for the two years is 
apparent. 

FigurelHlshows that our experiment was stable. We choose 
to interpret the excess as due to gamma-rays from this 
well-known emitter. We use the same acceptances as for 



Mrk 421, again matched to the conditions of each data 
run. The flux at 110 GeV thus obtained for the year 2000 is 
(0.5 ±0.15) X lO"!" ergcm-^s-i (Statistical errors only). 
Combining both years decreases the result by less than 
la. Figure 1131 shows the flux superimposed on an SED 
taken from ( Kataoka et~ar 1999^1 . The ASCA X-ray and 
EGRET gamma ray data were simultaneous with each 
other. We have added the BeppoSAX X-ray data taken at 
the same time as CELESTE'S. The TeV spectrum for the 
1998 season from the CAT Cherenkov imager corresponds 
to a quiet X-ray period l|Piron et al. 200 lb| . with an aver- 
age of 0.6 ASM counts per second, close to the year 2000 
rate, and no significant flares. 

6.3. lES 1426+428 

Table and the light curve in Figure El show that 
CELESTE detected no gamma-rays from lES 1426±428, 
indicating that the 2a upper flux limit at 80 GeV is 
0.2 X 10-1° ergcm^^s^i. The average ASM X-ray flux 
during CELESTE'S observations was 0.2 ±0.1 counts per 
second. For the two seasons leading to the detection of this 
blazar with the Whipple telescope the ASM rate ranged 
from 0.15 ± 0.04 to 0.65 ± 0.06, with an average near 0.3, 
with no discernible correlation between the weak X-ray 
and TeV signals l|Horan et al. 2002|l . 

Prior to the recent HESS result in 
(|Aharonian et al. 20061) . lES 1426±428 was the most 
distant blazar used to constrain the density of extragalac- 
tic diffuse infrared light, possible because TeV photons 
can be absorbed by low energy photons via electron- 
positron pair production, 77 — > e^e~. In particular. 



10 



CELESTE: Mrk 421, Mrk 501, and lES 1426+428 at 100 GeV 



51960 



51970 



51980 



51990 



52000 



52010 




30 



20 



- 10 



-10 



- 15 



10 



- 5 



E 
"S 
E 
E 
ns 
oi 



E 

E 
E 

(0 
O) 



51950 



51960 



51970 



51980 



51990 



52000 



52010 



52020, 



(MJD) 



Fig. 8. Mrk 421 daily averages in 2001. Epoch (2) is defined in Table El and discussed in the text. CELESTE data 
runs fro m MJD 51957 to 52015 (February 17 to April 16). STA CEE data is from (B oone et al. 200 2). HEGRA data 
is from (jAharonian et al. 2002|l . and the VERITAS data is from (|Holder et al. 200T)l . The other data points are as in 
the preceding figure. Note the different left- and right-hand axes. 



Data set 



Cut 



Number of events 

iVon iVoff A^On - 



A^off 



Significance 



Year 2000 (8.3 h) 



Raw data 
Software trigger 
All cuts 



616 194 
362 863 
73 034 



613 337 2 857 
359 914 2 949 
71518 1516 



Year 2001 (2.0 h) 



Raw data 
Software trigger 
All cuts 



130 820 
67040 
13 957 



131 584 
67410 
14 059 



-764 
-370 
-102 



All data (10.3 h) Raw data 747 014 744 921 2 093 

Software trigger 429 903 427 323 2 580 
All cuts 86 991 85 578 1413 

Table 4. Number of events before and after cuts for the Mrk 501 data sets. 
~ 80% data acquisition efficiency. 



2.9 
3.4 



-0.9 
-0.6 



2.4 
2.9 



Excess 
(min-i) 



5.9 
3.1 



-3.1 
-0.85 



4.2 
2.3 



Nc is calculated after correcting for the 



IjDwek fc Krennrich 2005( ) took the spectra observed at 
high energies and, for a range of plausible assumptions 
about infrared spectral shapes and densities, removed 
absorption effects. Extrapolations of the de-absorbed 
spectra to 100 GeV suggested fluxes detectable by 
CELESTE. In the end, our non-detection during a 



low X-ray state does not constrain the diffuse infrared 
densities. 

7. Conclusions 

We have provided some of the first astrophysical flux mea- 
surements around 100 GeV, beyond the reach of EGRET 



CELESTE: Mrk 421, Mrk 501, and lES 1426+428 at 100 GeV 



11 



53050 53055 53060 53065 



53070 53075 



53080 



E 
E 
<o 

O) 
O 



1 - 




8.0 



E 

i 4.0 



2.0 



0.0 



iO Asr 


1 1 1 1 

\/\ (2-12 


ke 


1 1 1 1 

V) 


1 1 1 1 


1 1 1 1 


1 1 1 1 


1 1 1 1 


1 ' ' 

▲ F 


CA: 


















<3-20 


- i- 








A 








*■ * 














----XI] 


D 

m 


t 

_L 


^ - 
-1— 1-^ 


" 1 1 1 1 - 




c 




;o CE 


III 
.ESTE 


-1- 


.lOOGeV) 




1 1 1 1 


1 1 1 1 


Mil- 


|— 1— 1- 




- 


® 














- ^ 


- 

- 
- 


■ Mil- 

i-e VEF 


RITAS 


E> 










Mil- 

^- 




-1— 1-: 


300GeV) 










1 


-j 
























m 




----^ 
















3 

1 . , , 


c 

1 


J® 
, , , , 


, , , , 


, , , , 






1 . . 





60 
40 
20 




53045 53050 



53055 



53060 53065 



53070 53075 



53080 53085, 



(MJD) 



Fig. 9. Mrk 421 daily averages in February and March 2004. Epochs (4) and (5) are defined in Table 3 and discussed 
in the text. The RXTE PCA and VERITAS data are taken from ( |Blazejowski et al. 2005) 1, while the ASM data comes 
from Ijsee ASM in references|l . 



Number of events 



Data set 



Cut 



Non 



N, 



Non - N, 



Significance 



All data (SPV, 8.8 h) Raw data 675 340 673 239 2 101 

Software trigger 455 794 453 660 2 134 2.0 
All cuts 64 230 64 433 -203 -0.5 

Table 5. Number of events before and after cuts for the lES 1426+428 data sets. N^^ is calculated after correcting 
for the ~ 80% data acquisition efficiency. 



on the Compton CGRO and below the range of previ- 
ous ground-based telescopes, by recording atmospheric 
Cherenkov light with the reconverted solar tower facil- 
ity at Themis. We have increased our detector sensitiv- 
ity as compared to our earlier work, to a level near that 
predicted in our original proposal. The gain came mainly 
from a better analysis procedure but also from refinements 
in our detector simulation and electronics improvements. 
80 % of the gamma-rays we detect have energies between 
50 and 350 GeV. Bad weather at our site lead to a low 
duty-cycle, resulting in a modest number of detections in 
spite of our detector's good performance. 



Our results for the blazar Mrk 421 help constrain the 
shape of the high energy peak of the spectral energy dis- 
tribution. We have made the first detection of Mrk 501 in 
this range, although the detection is weak, and we have 
placed an upper limit on the flux from the distant blazar 
lES 1426+428 during a period of X-ray inactivity. 
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Fig. 10. CELESTE average flux measurement (gray bowtie) for Mrk 421 for 2001 February (MJD 51956.98 to 51963.11, 
Epoch (1) in Table 3 and text), superimposed on the spectral energy distribution taken from the multiwavelength 
study by (IBlazejows ki et al. 2005| ). The envelope of the bowtie corresponds to the assumptions of a = 1.8 and 2.2 
differential spectral indices. The inner set of lines above and below the bowtie are the statistical errors, while the outer 
set has in addition a 25% systematic uncertainty added in quadrature. The ASM X-ray activity level in February 
2001 corresponds to the "medium" activity state (the lower of the solid-curved high energy peaks) and thus to the 
VERITAS TeV measurements represented by the squares. The CAT spectrum in ( Piron et al. 2001a() for the year 2000 
(thick black line) gives the same integral flux as for the nights in February 2001. 
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Fig. 12. Gamma- and X-ray light curves for Mrk 501. MJD 51600 and 52000 correspond to 2000 February 26 and 
2001 April 4. TOP: ASM 5-day running averages of tlie 2 to 12 keV X-ray flux, in cts/s Ijsee ASM in references|l . 
and BeppoSax 2 to 12 keV fluxes in IQ-^*' ergcm"^ g-i, from (jMassaro et al. 2004|l . BOTTOM: CELESTE daily flux 
averages. 
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Fig. 13. CELESTE flux measurement for Mrk501 from tlie year 2000 (bowtie), bounded by the assumptions of a = 1.8 
and 2.2 differential spectral indices. The inner set of lines, above and below the bowtie, correspond to the statistical 
uncertainty, while the outer set have in addition the 25% systematic uncertainty added in quadrature. The Beppo 
X-ray data is nearly simultaneous IjJVlassaro et al. 2004|l . The EGRET data and the ASCA X-ray measurements are 
simultaneous with each other, and correspond to the SSC model curves IjKataoka et al. 1999|l . The TeV spectra are 
from the CAT Cherenkov imager HPiron et al. 2001b)l . 
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Fig. 14. Light curves for lES 1426+428 - CELESTE detected no signal. The CELESTE daily averages run from 2002 
March 22 to 2004 April 25. The ASM points are 5 day running averages. 



